i = 0
nn.results <- list()
size_seq <- seq(4,12,2)
# Loop over years
for (year in start.year:end.year) {
  i = i + 1
  print(year)
  # Generate dependent variable
  Violence.mean <-  aggregate(x = as.matrix(dta[dta[,t.var]< year,
                                                dv]),
                              by = list(dta[dta[,t.var]< year,
                                            loc.var]),
                              FUN = mean)
  colnames(Violence.mean) <- c(loc.var,
                               "dv_mean")
  Violence.new <- merge(as.matrix(dta[dta[,t.var]< year,
                                      c(loc.var,t.var,dv)]),
                        Violence.mean,
                        by = loc.var)
  Violence.past <- Violence.new[,dv]-Violence.new$dv_mean
  # Get optimal parameters from cross validation
  # Fit Neural Network
  nn.results[[i]] <- list()
  dta.past <- dta[dta[,t.var]< year,]
  N <- length(dta.past[,dv])
  cv_results <- foreach (cv=1:cv.runs, 
                         .combine='rbind',
                         .packages = c("nnet")) %dopar% {   # Loop over CV runs
                           set.seed(52*cv) 
                           shuffle <- sample(1:N,
                                             N,
                                             replace=FALSE)
                           predictions <- matrix(NA,nrow=N,ncol=length(size_seq))
                           sampleSplit = c(0,round((N/5*c(1:5))))
                           
                           for (f in 1:5) {
                             # Remove constant variables from predictors
                             train.RHS <- dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                        rhs]
                             test.RHS <- dta.past[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                       rhs]
                             # Remove variables that are constant
                             not.constant <- apply(train.RHS,2,sd)!=0
                             train.RHS <- train.RHS[,not.constant]
                             test.RHS <- test.RHS[,not.constant]
                             # Find PC weights
                             pcs.fit <- prcomp(train.RHS,center = TRUE, scale = TRUE)
                             train.pcs <- predict(pcs.fit,
                                                  train.RHS)[,1:min(nn_demean_params$pcNum,
                                                                    ncol(train.RHS))]
                             test.pcs <- predict(pcs.fit,
                                                 test.RHS)[,1:min(nn_demean_params$pcNum,
                                                                  ncol(train.RHS))]
                             # Fit Neural Network
                             for (s in 1:length(size_seq)) {
                               # Fit Neural Network
                               cv.fit = nnet(x = train.pcs,
                                             y = as.matrix(Violence.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])]]),
                                             size = size_seq[s], 
                                             decay = nn_demean_params$decay,
                                             trace = nn_demean_params$trace,
                                             linout = nn_demean_params$linout
                               )
                               predictions[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                         s] <- predict(cv.fit,
                                                       newdata = as.matrix(test.pcs))
                             }
                           }
                           cv_results <- c("CV Run"=0,
                                           "size"=0,
                                           "MSE"=100)
                           for (s in 1:length(size_seq)) {
                             MSE <- mean((Violence.past - predictions[,s])^2)
                             if (MSE<cv_results["MSE"]) {
                               cv_results <- c("CV Run"=cv,
                                               "Size" = size_seq[s],
                                               "MSE" = MSE)
                               best.s <- s
                             }
                           }
                           return(cv_results)
                         } 
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  # Fit Neural Networks
  nn.results[[i]]$size <- round(median(cv_results[,"Size"]))
  dta.past.rhs <- dta.past[,rhs]
  
  not.constant <- apply(dta.past.rhs,2,sd)!=0
  dta.past.rhs <- dta.past.rhs[,not.constant]
  pcs.fit <- prcomp(dta.past.rhs,
                    center = TRUE, 
                    scale = TRUE)
  insample.pcs <- predict(pcs.fit,dta.past.rhs)[,1:min(nn_demean_params$pcNum,ncol(dta.past.rhs))]
  outsample.pcs <- predict(pcs.fit,dta[dta[,t.var]== year,])[,1:min(nn_demean_params$pcNum,ncol(dta.past.rhs))]
  
  set.seed(i)
  nn.results[[i]]$mod <- nnet(x = insample.pcs,
                              y = Violence.past,
                              size = nn.results[[i]]$size,
                              decay = nn_demean_params$decay,
                              trace = nn_demean_params$trace,
                              linout = nn_demean_params$linout
  )
  nn.results[[i]]$fit.oos <- predict(nn.results[[i]]$mod,
                                     newdata = outsample.pcs)
  Violence.new <- merge(as.matrix(dta[dta[,t.var]== year,
                                           c(loc.var,t.var,dv)]),
                        Violence.mean,
                        by = loc.var)
  nn.results[[i]]$actual.oos <- Violence.new[,dv]-Violence.new$dv_mean
  print(year)
}
save(nn.results,file=paste(modeldir,"/",
                            country,
                            "_nn_",
                            "demean",
                            "_",
                            rhs.group,fileext,
                            ".RData",
                            sep = ""))